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MATRIX-EQUATION-BASED STRATEGIES EOR 
CONVECTION-DIFFUSION EQUATIONS* 

DAVIDE PALITTAt AND VALERIA SIMONCINit 


Abstract. We are interested in the numerical solution of nonsymmetric linear systems arising 
from the discretization of convection-diffusion partial differential equations with separable coefficients 
and dominant convection. Preconditioners based on the matrix equation formulation of the problem 
are proposed, which naturally approximate the original discretized problem. For certain types of 
convection coefficients, we show that the explicit solution of the matrix equation can effectively 
replace the linear system solution. Numerical experiments with data stemming from two and three 
dimensional problems are reported, illustrating the potential of the proposed methodology. 
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1. Introduction. We are interested in the numerical solution of the convection- 
diffusion partial differential equation 


eAu -I- w • Vm = /, in n C 


( 1 . 1 ) 


with d = 2,3, where w is the convection vector, while e is the viscosity parameter. 
In particular, we consider the dominant convection case, that is |w| ^ e, and we 
assume that w is incompressible, that is div(w) = 0. Moreover, we assume that the 
components of w are separable functions in the space variables. For simplicity the 
equation is equipped with Dirichlet boundary conditions; the analyzed procedures 
could be used with Neumann boundary conditions as well. 

Standard finite difference or finite element discretizations yield the algebraic large 
nonsymmetric linear system 


Au = f, with A e 


( 1 . 2 ) 


The discretization phase is crucial to obtain a reliable numerical solution to (1.1), as 
spurious oscillations may occur in the approximate solution in the dominant convec¬ 
tion regime; see, e.g., [10]. A variety of strategies is available, which either take a fine 
enough discretization, or appropriately modify the differential operator so as to limit 
the convection side-effects [14]; see [30] for a recent essay on the subject. In our nu¬ 
merical experiments we used the former approach, though the proposed methodology 
could also be adapted to handle, e.g., SUPG techniques [23]. 

We are interested in exploring preconditioning strategies for solving (1.2), that 
stem from knowledge of the differential operator and its discretization. In particular, 
a computationally cheap approximation to the original differential operator will be 
employed, that exploits the original matrix equation structure of the problem. The 
idea of using a simplified operator as a preconditioner is well known in the convection- 
diffusion equation literature. More precisely, classical strategies use the symmetric 
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part of the operator - see, e.g., [13],[1], and in particular the discussion in [19] - due 
to the fact that a symmetric problem is usually cheaper to solve than the original 
nonsymmetric system if fast solvers such as fast Fourier transform, cyclic reduction 
or multigrid can be employed. In the case of nonsymmetric preconditioners, attempts 
have focussed on simplifying the operator so as to still use fast solvers, and this is 
achieved, e.g., by imposing constant coefficients or by dropping some of the first order 
terms; see, e.g., [2],[7],[9],[13]. However, their effectiveness and computational cost 
of their application have not been clearly assessed. Finally, most publications aim 
at analyzing the two-dimensional (2D) problem, whereas the three-dimensional (3D) 
problem represents a challenging task, especially from a computational point of view. 

We exploit the matrix structure of the discretized problem to construct a cost- 
effective nonsymmetric approximation: the application of the new preconditioner 
amounts to solving a Sylvester matrix equation, whose coefficient matrices are tightly 
related to the original discretized problem. Other authors have used matrix equations 
either as a solver or as a preconditioner for the system stemming from (1.1); see, e.g., 
[32], [33], [29], [26]. However, on the one hand, only very simplified models have been 
considered, on the other hand, to the best of our knowledge no performance evalu¬ 
ation has ever been performed with respect to more common multilevel techniques. 
We derive a new preconditioning operator by first writing down the full multiterm 
matrix equation corresponding to finite difference discretization of the two or three- 
dimensional problem and then by appropriately simplifying the operator; then we 
iteratively solve the associated Sylvester equation so as to achieve good performance 
of the overall preconditioning phase. The choice of the Sylvester equation solver is 
crucial for assessing the preconditioning costs. In the literature, ADI and cheap solvers 
for the Kronecker formulation were suggested for various (e.g. symmetric) variants of 
the problem [7],[15],[33]; here we exploit a recently developed efficient iterative solver 
for the Sylvester equation, whose cost per inner iteration for a 2D problem may be 
as low as 0(n), where n is the one-dimensional problem size; this solver allows us 
to preserve the leading coefficients of the two first order derivatives without losing 
efficiency. 

Our preliminary numerical results are promising and clearly illustrate the poten¬ 
tial of this approach: the overall solution compares rather well with state-of-the-art 
and finely tuned algebraic multigrid preconditioners. When used as a solver, in both 
2D and 3D problems, its performance is superior to other available approaches. Fi¬ 
nally, our numerical experiments confirm theoretical results in the literature on the 
robustness of the approach with respect to the mesh parameter, namely the number 
of iterations of the system solver is essentially independent of the problem size. 

Our aim is to derive a proof-of-concept preconditioning procedure for the finite 
difference discretization of the problem on rectangular and parallelepipedal domains, 
as a first step towards its higher level development for the finite element discretization 
on more general domains. 

Here is a synopsis of the paper. In section 2 we derive the matrix formulation of the 
two-dimensional discretized problem, from which the standard form (1.2) is derived 
by means of Kronecker product expansion. Section 3 describes how the boundary 
conditions can be embedded in the matrix formulation. Section 4 uses the previously 
derived form to define the preconditioning operator and its associated Sylvester matrix 
equation, while implementation details are given in section 4.1. The procedure is then 
generalized to three dimensional problems in section 5, while numerical experiments 
with systems stemming from both two and three dimensional data are reported in 
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section 6. Finally, our conclusions are given in section 7. 

2. A matrix oriented formulation. In this section we reformulate the alge¬ 
braic problem (1.1) in terms of a multiterm linear matrix equation. This derivation 
will be used to introduce our preconditioner, to be applied to a Krylov subspace 
method as an acceleration strategy for solving (1.2). For the ease of presentation, we 
shall first concentrate on the two-dimensional problem, and then extend our derivation 
to the three-dimensional case in section 5. 

We start by recalling the matrix equation associated with the discretization by 
five-point stencil finite differences of the Poisson equation —Au = / on a rectangular 
domain fl C For the sake of simplicity, we shall assume that il = (0,1)^. Let fl/i 
be a uniform discretization of fl, with nodes {xi, yj), i,j = 1,..., n — 1. Then assum¬ 
ing homogeneous Dirichlet boundary conditions are used, centered finite difference 
discretization leads to the linear system (1.2) with 


A = r„_i 0 /„_i -b In-i 0 r. 


and Tn-i = tridiag(—1,2, —1) g is the symmetric tridiagonal matrix 

approximating the second-order derivative in one-dimension, while the entries of u 
contain an approximation to u at the nodes, having used a lexicographic order of the 
entries. 

We thus take a step back, and describe in details the process leading to the Kro- 
necker formulation, with the aim of deriving its matrix counterpart. This description 
will allow us to also include the boundary conditions in a systematic manner. 

Let VLfi be a uniform discretization of the closed domain fl, with equidistant points 
in each direction, {xi,yj), i, j = 0,... ,n. Analogously, Uij = U(xi,yj) is the value of 
the approximation U to it at the nodes. For each i,j = 1,..., n — 1 we have the usual 
approximations 


^XX (^2 ; yj ) 



and analogously for the y direction, but from the right, 


Uyy 




Uij-i — 2Uij -b Uij^i 



Let 


/ * 


* * 


\ 


* -2 1 



•. 1 

1 -2 = 1 = 


g R("+l)x(n-Hl), (^2.1) 


V 


* / 


the unspecified values are associated with boundary values of U and will be 
discussed in section 3. Collecting these relations for all rows I’s and for all columns 
j’s, for the whole domain we obtain 
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With these approximations we can write the following classical matrix form of the 
finite difference discretization of the Poisson equation on a square domain (see, e.g., 

[31]) 

TU + UT = F, where Fij = f{xi,yj) + b.c.. (2-2) 

Except for the boundary conditions, the Kronecker formulation of (2.2) gives the same 
form as (1.2). 

For the convection-diffusion equation with separable coefRcients a similar deriva¬ 
tion provides a multiterm linear matrix equation. We state the result in the following 
proposition, where separable convection coefficients are assumed. To this end, we 
define the matrix 



* 

0 

-1 


1 


-1 


1 

0 

* 


g ]g("’+l)x("+l). 


* 

* 


(2.3) 


which represents the centered finite difference approximation of the first order one 
dimensional (ID) derivative on a uniformly discretized interval. 

Proposition 2.1. Assume that the convection vector w = {wi,W 2 ) satisfies 
wi = 4>i{x)'ilJi{y) and W 2 = (t> 2 {x)ii 2 {y) ■ Let {xi,yj) € Clh, ij = 0,... ,n and set = 
diag((()fc(a;o),..., ((>fc(a;„)) and 4'^ = diiig{ipkiyo), ■ ■ ■ ,i’kiyn)), k = 1,2. Then with 
the previous notation, the centered finite-difference discretization of the differential 
operator in (1.1) leads to the following operator: 

Ch : U ^ eTU + eUT + (4>iB)t/4'i -h $2^(5'^ 4 - 2 ). (2.4) 


Proof. The first two terms of Ch{U) were derived for (2.2). We are left with 
showing that the first order term can be expressed in terms of the ID discretization 
matrix B in (2.3). We have 


(Pi[Xi)tpi{yj)u,,{xi,yj) ze cpffxi) - — - Fiiyj) 




Ui — lj 
Ui+ij 


V’i(yj), 


and analogously, 


(l)2iXi)lf'2{yj)Uy{Xi,yj) 


2j^4’2{xi)[Uij—i, Uij, Uij+i] 


-1 

0 

1 


^2(yi)- 


Collecting these results for all grid nodes and recalling that Uij = U{xi,yj), we obtain 


{<f2i.Xi)llj2{y3)Uy{Xi,yj))^ .^^^ ^^ ^ 


$2DS^4'2, 
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and the result follows. □ 

The operator in (2.4) is a linear multiterm matrix function, and it is a general 
variant of the matrix equations already present in the early literature on difference 
methods [4]. This type of equations now often arises in the discretization of partial 
differential equations with stochastic terms; see [28] for a description of this and 
other contexts where these operators arise. In the deterministic setting, however, 
the difficulty of explicitly dealing with all terms have led the scientific community to 
abandon this form, in favor of the Kronecker formulation. 


3. Imposing the boundary conditions. The algebraic problem needs to be 
completed by imposing the boundary values. These will fill up the undefined entries 
in the coefficient matrices, and in the right-hand side matrix. To this end, we recall 
that with the given ordering of the elements in U, the first and last columns, Uei and 
Ucn+i resp., correspond to the boundary sides y = 0 and y = 1, whereas the hrst 
and last rows take up the values at the boundary sides a; = 0 and a; = 1. With this 
notation, we wish to complete the corners of the matrices T and B, giving rise to the 
matrices Ti, T 2 and Bi, B 2 , respectively, so that the following matrix equation is well 
defined for (xi,yj) € Cth'- 

eTiU + eUT2 + $1^1 C/4'1 -f 4>2C/B2^2 = F. (3.1) 


With the same notation as for [/, the entries of F corresponding to i,j € {0,n} will 
contain contributions from the boundary values of C/, which are determined next. 

For the boundary conditions to be satisfied, the operator Ch{U) = eTiU + eUT2 + 
<I>ii?iC/4'i -I- ^2UB2'i'2 should act as the (scaled) identity operator for points at the 
boundary. To this end, from the generic matrix T we define the matrix Ti as follows: 


Ti = - 




/ -1 0 

1 -2 1 

1 -2 1 


V 


1 


(n-|-l)x(ri+l) 


0 -1 / 

while the matrix corresponding to the first order operator (B in the generic case) can 
be written as 

/ 0 0 \ 

-1 0 1 
-1 0 1 


2h 


-1 0 1 

0 0 / 


(ri+l)x(n-|-l) 


For the derivative in the y direction, right multiplication should also act like the 
identity, therefore we can define T 2 = and correspondingly, B 2 = Bf. We are 
thus ready to define the missing entries in F so that (3.1) holds. For the first column, 
that is for the side y = 0, we write 


Fei = {fTxU -f CUT2 ^\B\U^\ $21^1^2^2)61 

= eT\U61 -h eC/T2Ci -f ^iBit/^liei -I- $2!! 132 $261 
= eTiC/ei -I- -^Uei + ^i(yo)^iBiUei, 
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where we used the fact that = 0- Similar reasonings ensure that the boundary 

values at y = 1 are imposed, thus defining Fe„+i. For the side a: = 0 we have 

elF = e^ (eTit/ + ef/Tz + $iSi C/4'1 + ^2UB2'l>2) = e{F 
= CieTiU + efeUT2 + ef<I>ii?iC/4'i + ei^2UB2'^2 

= -^elU+ ee{UT2 + c^2{xo)elUB2'^2, 

where the fact that = 0 was used. The definition of e^^iF follows analo¬ 

gously. 

We stress that maintaining the same boundary conditions in the designed precon¬ 
ditioner has other benefits, in addition to accurately reproducing the original operator. 
Indeed, it was shown in [20] that if the preconditioner is itself a matrix arising from 
the discretization of an elliptic operator, then under certain regularity assumptions 
the L 2 -norm of the preconditioned matrix AP~^ is uniformly bounded with respect 
to the discretization parameter. 

4. The new preconditioner. The main obstacle in directly solving equation 
(2.4) is given by its many general terms. Indeed, had the equation terms with co¬ 
efficient matrices on either side of the unknown matrix, the problem would be read¬ 
ily recognized as a matrix Sylvester equation. This scenario does occur whenever 
some of the convection terms have a simplified structure, so that those terms only 
depend on some of the variables. This is the case, for instance, with the opera¬ 
tor jC{u) = —Am -I- (j)i{x)ux, which gives rise to Ch{U) = TiU + UT 2 + ^iBiU = 
(Ti -I- ^iBi)U + UT 2 . In this case, computational strategies based on the matrix 
equation can be employed; see, e.g., [32],[29], and [8] and its references for some early 
methods; see also section 6 for some examples with state-of-the-art approaches. In 
the general multiterm case, a simplified version of the left-hand side in (2.4) can be 
used as an acceleration operator, in the form of a matrix equation solver; a similar 
strategy was adopted in [7], with encouraging numerical results with early solution 
methods. Let 


£h-Y^ eTiY + eYT2 + ($iBi)y^i + <^>2Y{B2'S2) 

be the operator associated with (3.1). Then Ch can be approximated by replacing 
two of the diagonal coefficient matrices, namely 

^ 2 ^ hi, (4.1) 

where the scalars ipi and <j )2 represent some average of the functions and (f )2 on the 
given domain. In our experiments we used 

i i 

but other strategies may be considered. Substituting the approximations (4.1) in Ch, 
the matrix Y can be collected, yielding the following approximation to Ch- 

V:Y^{eT + iPi^iB)Y + Y{eT + (4.2) 

The use of V in an acceleration context consists of applying to a given matrix 
G. This corresponds to solving the following Sylvester equation 


V{Y) = G O {eT^Fh^^B^)Y+ Y{eT2 + B2'i!2h) = G. (4.3) 




matrix-equation strategies for convection-diffusion equations 


7 


This matrix equation has a unique solution for any G 7 ^ 0 if and only if the spectra of 
eTi+'0i$i_Bi and — (eT 2 +have no common eigenvalues. This is ensured for 
instance if both coefficient matrices in (4.3) have eigenvalues in C“'". The numerical 
solution of (4.3) will be discussed in section 4.1. 

We would like to point out that the approximation in (4.1) corresponds to iirtro- 
ducing the following modified convection vector 

w := (r^i(/)i(a;),^ 2 ^/’ 2 ( 2 /)) ~ w. (4.4) 

so that the continuous operator C{u) = —eAu + w ■ Vm is approximated by C{u) = 
—eAu + w ■ Vu. Similar preconditioning strategies relying on simplified (non-self- 
adjoint) versions of the operator C have been developed in the past, see, e.g., [ 2 ], 
[7], [9], [13]; some of them use a matrix equation oriented formnlation, which allows 
one to significantly lower the computational cost, compared with the Kronecker form, 
while keeping non-constant coefficients. In all these cases, however, no performance 
assessment has been clearly analyzed, and in fact, the overall cost heavily depends 
on the complexity of the method used to apply the preconditioner. By employing a 
matrix equation formulation, the computational cost is in general significantly lower 
than if one were to solve at each iteration a nonsymmetric linear system that is very 
close to the original one (unless special strategies can be designed for the latter), see, 
e.g., the discussion in [2, page 1514]. We refer to [23] for experimental comparisons 
with inner-outer procedures that explore incomplete LU preconditioners. 

In the continuous case, our strategy determines a preconditioner that is^ “com¬ 
pact - equivalent” to the original operator C [2, Proposition 3.5]. In particular, this 
property implies that if one were to use the conjugate gradient method to the normal 
equation associated with our preconditioned problem, its convergence rate would be 
bounded by a quantity that is independent of the mesh parameter [2, Theorem 4.7]. 
In our numerical experiments we did not employ the normal equation, but we can still 
empirically infer that mesh independence is maintained. Finally, we should mention 
that this property does not necessarily imply that the preconditioner is good, that is 
that a low number of iterations is obtained. 

4.1. Implementation details. Since A in (1.2) is nonsymmetric, the Krylov 
subspace solver GMRES ([25]) for nonsymmetric problems can be used, and the oper¬ 
ator preconditioner 'P~^ is applied from the right. In a standard implementation, at 
each iteration k of GMRES the preconditioner is applied to a vector as yu = V~^gk, 
where gk is a vector with N = (n + 1)^ components. In our matrix strategy, we first 
transform the vector gk into the matrix Gk such that gk = vec(Gfc), and then obtain 
yk = vec{Yk) as the solution to the Sylvester equation 

(eTi -|- 'tjji^iBi)Y + Y{eT2 + 52^2^' 2 ) = Gk, (4-5) 

where the coefficient matrices have both dimensions n-l-1. For n up to a few hundreds, 
this equation may be efficiently solved by means of the Bartels-Stewart method [3], 
whose computational cost is G((n-|-1)^). This cost should be compared with the cost 
of solving an inner system of size (n-|-1)^ x (n-l-1)^ by either a sparse direct solver, or 
by a fast solver if applicable. If a different number of nodes is used in each direction, 
or a three-dimensional problem is solved, then the two coefficient matrices will have 

^Roughly speaking, two elliptic operators are compact - equivalent if their principal parts coincide 
up to a constant factor, and they have homogeneous Dirichlet conditions on the same part of the 
boundary. 
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different size, without any significant impact on the solution process, though the cost 
will change accordingly. 

For larger n, the application of the Bartels-Stewart method may become too 
expensive to be competitive, therefore an iterative method should be employed. In 
this case, the Sylvester equation is solved up to a certain tolerance, so that only 
an approximation to V is obtained. To make the strategy cost and memory effective 
when using state-of-the-art iterative methods, the right-hand side should be low rank. 
Therefore, Gk is also approximated by a low-rank matrix by means of a truncated 
singular value decomposition. The overall effectiveness of the acceleration process 
thus depends on how much information the truncation retains, and how accurate the 
iterative solution will be. These require the tuning of three parameters (the truncation 
and stopping tolerances, and the maximum rank allowed for Gk) that we set a-priori in 
all our numerical experiments: tol_truncation=10“^, tol_inner=10“'^, r_max=10. 
The performance did not seem to be influenced by small variations of the values of 
these parameters. 

Finally, we mention that if an iterative solver is used for (4.3), then the precondi¬ 
tioning step makes the overall process a nonlinear operation, as the inexact operator 
“changes” at each system solver iteration. As a consequence, the flexible ver¬ 
sion of GMRES, named FGMRES in the following [24], is used, giving rise to an 
inner-outer iteration. 

The iterative solution of the Sylvester equation is performed by means of the 
Extended Krylov subspace method (KPIK), originally proposed in [27] for the Lya¬ 
punov equation, and then adapted to solving the Sylvester equation in [6]. The relative 
residual norm is used as stopping criterion, and the residual is checked at every other 
iteration, as suggested in [27]. Other efficient Sylvester equation solvers could be used, 
see the recent survey [28]. However, in a preconditioning framework, the fact that 
the coefficient matrices can be factorized once for all makes KPIK very appealing in 
terms of computational costs; see [27] for further details. 


5. The three-dimensional case. The 3D convection-diffusion equation can be 
stated as in (1.1), for H C To convey our idea, we again first focus on the 
Poisson equation, and then generalize the matrix formulation to the hnite difference 
discretization of the non-self-adjoint problem (1.1). For the sake of simplicity, we shall 
assume that D = (0,1)^, though more general parallelepipedal domains could also be 
considered. We discretize H with equidistant nodes in each direction, (xi,yj, Zk), 
for i,j, fc = 0,... ,n. To fix the ideas, let = U{xi,yj,Zk) denote the value of the 
approximation U to it at the node (xi,yj,Zk) (other orderings may be more convenient 
depending on the equation properties). We also define the tall matrix 


t/(o) 


n 


u = 


[/(") 



where Cj is the jth canonical vector of 
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Let T be as defined in (2.1). Then, for I G the identity matrix, 

n n 

-U,. - ^(efc+1 0 T[/«) = (I0 T) ^(efc+i 0 C/^) = (/ 0 Tp, 

k^O k^O 

n n 

-uyy ~ Y.^ek+1 0 U^'^'^T) = ® 

k^O k^O 

-Uzz ~ {T <^I)U. 

With these approximations we can thus obtain the following matrix form of the 
finite difference discretization of the Poisson equation: 

il0T)U+UT + {T0l)U = F, (5.1) 

where F = J2k=oi^k+i®F^^'>) G and iF^'^'>)ij = f{xi,yj,Zk). The Kro- 

necker formulation of the matrix equation (5.1) determines the usual approximation 
of the Laplacian operator by seven-point stencil finite differences, 

Aw/0/0T-b/0T0/-bT070/G 

For the convection-diffusion equation with separable coefficients a similar derivation 
provides a multiterm linear matrix equation. We state the result in the following 
proposition. 

Proposition 5.1. Assume that the convection vector w = {wi,W 2 ,w^) sat¬ 
isfies Wi = cj)i{x)fi{y)vi{z), W 2 = 'p 2 {x)i> 2 iy)v 2 {z), and W 3 = 4>3{x)'ip3{y)v3{z). 
Let {xi,yj, Zk), i,j,k = 0,...,n be the grid nodes discretizing Q with mesh size h, 
and set = dia.g{<pi{xo),..., (j)i{xn)), = diag(V'^(yo), ■ • •, i/’^(2/n)), and Ti = 

diag(p^(2;o), ■ • •, vd^n)), ^=1,2,3. 

Then, with B as defined in (2.3), the centered finite-difference discretization of 
the differential operator in (1.1) leads to the following operator: 

Ch ■■ (/0T)W-fWr+(r®/)W+(Ti®$iB)W^i-H(T20$2)WB^^2 + [(T3B)04>3]Wvl'3. 

(5.2) 


Proof. The second order terms of Ch{U) correspond to a multiple of (5.1). We 
are thus left with showing that the first order term can be expressed by means of the 
ID discretization matrix B. We first fix A: = fc and we have 


(l)i{xi)'ilJi{yj)vi{zd)udxi,yj, zi) 


Analogously, 


U - u[xi-i,yj,zd . 

x’lyZkjfiii.Xi) 'ipi[yj) 


2h 




2h 

^1-1,3 

hJ 


diivj)- 


Mxi)'f2{yj)v2{zf,)uy{x^,yj, zd 


2h 


V2iz-,)Mx^)[US-lM^Ml 


.H-iJ 


-1 

0 

1 


d2{yj)- 


Collecting these results for all grid nodes {xi,yj, zd)i,j=o,...,n and recalling that = 
U{xi,yj,Zk), we obtain 

{(l>l{Xi)lpl{yj)vi{zj.)udXi,yj,zd)i,j=0,...,n - 'Ul(Zfc)^'lFC/(''%i, 







10 


Davide Palitta and Valeria Simoncini 


and 

{(i)2{Xi)lp2iyj)v2{Zk)Uy{Xi,yj,Zj,))ij=o^,,,^ri - V2{Zk)‘^2U^^^ 'i>2. 

Therefore, for all z nodes, 

{(j)i {Xt)ijji{yj)vi{zk)uxixi,yj,Zk))i,j,k=o,...,n 

n n 

^ [Ti ® /] ^(efc+i 0 = [Ti ® /](/ 0 $iB)[^(efc+i 0 

fc =0 fe =0 

= [Ti 0 /](/ 0 = (Ti 0 

and 


i4>2ix^)'^p2{yJ)v2iZ:^)Uy{Xi,yj,Zk))i,j,k=0,...,n 

n n 

^ [T2 0 /] Y.^ek +1 0 = [T2 0 /](/ 0 ^ 2 )[J 2 {ek+i 0 U^'^'>)]B^'i /2 

k^Q k^O 

= [T2 0 /](/ 0 <i>2)l(B^^2 = (Tz 0 <^?2pB^'^2- 

On the other hand, for the z direction it holds 


(1)3 (a;z) V'3 {yj)v3 {Zk)uz {Xi ,yj,Zk) 


V3{Zk)(l)3{Xi) 


u{xi,yj,Zk+i) -u{xi,yj,Zk-i) 


—v:i{zk)(t)3{xi)[-l,Q,l] 


2 h 

j-j{k — l) 

^ i,j 

i,j 

(fc+1) 




^,3 


Myj)- 


Collecting this relation for all blocks, 

n n 

^ (TgS 0 I) ^[efc+i 0 ($ 3 C/('=)^' 3 )] = (T3B 0 /)(/ 0 $ 3 )E(efe+l 0 t/W )]^3 

/ c =0 fc =0 

= (T3B 0 /)(/ 0 $3)^^^'3 = [(T3S) 0 $3]Z^«'3. 

and the result follows. □ 

Imposing the boundary conditions completely determines the entries of T in all 
three instances, as well as the missing entries in B. Following the same steps as for 
the 2D case, the matrix equation (5.2) can be written as 

{{mTi) + {T^ ®I)) U + U T3 + (Ti(g)<I-iBi) U 'ti + (T2(gi<I>2) U ^3^2 + [(TgSj) ® $3] ^>^3 = 0 

highlighting the presence of hve distinct terms in the matrix equation. With this 
ordering of the variables, it holds that B^ = B 2 and T 3 = T 2 . Note that a tensorial 
formulation could also be obtained by further “unrolling” the Kronecker products in 
some of the terms. All three directions of the tensor would have dimension n + 1, 
thus being the natural generalization of the 2D problem. In future research we will 
explore the possibility of explicitly solving the fully tensorized equation (with three 
terms), possibly using recently developed strategies [16],[18]. 

The matrix equation above may take the form of a standard Sylvester equation 
depending on whether some of the coefhcients and V(, i = 1,2,3 vanish; see 
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Example 6.5. In the generic case, a preconditioning operator can still be derived, e.g., 
by averaging the values of some of the coefficients, the way it was done in the 2D case. 
The new preconditioner then consists of an approximation of the equation (5.2). For 
instance, by approximating 4'i,tl'3 as 

^'3 ss ■03/, 

and also using $2 ® Tf 2 « y/, the following operator is obtained, 

V -.V ^ ((/0Ti) + (T2^ 0/)-tV;i(Tl0$lBl)-h03[(T3Sj’)0#3]) V-fV(r3+xS3'I'2) 

whose application entails the solution of a Sylvester equation with coefficient matrices 
of size (n + 1)^ x (n -I- 1)^ and (n -I- 1) x (n -|- 1), respectively. 

6. Numerical experiments. In this section several numerical experiments are 
presented using both two and three dimensional problems. Performance with respect 
to the problem parameter ~ the viscosity - and the discretization parameter - the 
meshsize - are considered. 

In both dimensional settings we first consider the case when the matrix equation 
framework can be used as a solver for the original equation: this corresponds to 
problems where the first order term coefficients only depend on the same variable as 
the corresponding derivative. Comparisons with either sparse direct solvers or with 
iterative solvers are shown. Then we report on our experience when using the matrix 
equation strategy as a preconditioner for more general convection-diffusion problems 
with separable coefficients. We compare the performance of the new approach with 
that of state-of-the-art algebraic multigrid preconditioners; experimental comparisons 
with (less performing) ILU-type preconditioners can be found in [23]. More precisely, 
we shall consider both the algebraic multigrid preconditioner MI20 [5] with GMRES 
as a solver, and AGMG [21] with flexible GMRES as a solver (in AGMG all default 
parameters were used, while in MI20, control. one_pass_coarsen was set to 1). Both 
strategies have been shown to be applicable to convection-diffusion equations, and in 
particular, for AGMG it was shown in [22] that this variable preconditioning strategy 
is well suited for both 2D and 3D problems. 

We stress that our experimental comparisons somewhat penalize our approach, 
as the other preconditioners are in fact fortran90 fully compiled codes, for which a 
mex file was made available. Their performance is thus expected to be superior to 
interpreted Matlab functions, on which our preconditioner is based. Nonetheless, the 
reported results show that the new strategy is still competitive. 

Except for Example 6.3, all problem data were obtained by centered finite dif¬ 
ference discretization on the given domain; for all experiments the grid fineness was 
chosen so as to avoid spurious oscillations in the numerical solution for the coarsest 
grid used. 

All experiments were performed with Matlab Version 7.13.0.564 (R2011b) on a 
Dell Latitude laptop running ubuntu 14.04 with 4 GPUs at 2.10GHz. 

6.1. Two-dimensional problems. Whenever the convection vector has the 
simplified form w = (wi,W 2 ) = 1 ^ 2 ( 1 /)), that is ipi and 02 are constant func¬ 

tions, the matrix formulation (2.4) reduces to 

eTiU -f eUT2 + (4>iHi)t/ -k U{B2'i>2) = F, 

which is a Sylvester matrix equation, and can thus be solved directly, with no further 
approximation. In the first two examples, we thus consider equations leading to this 
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simplified form, and we compare the performance of either two direct solvers - in the 
matrix and vector equation regimes respectively, or of two iterative solvers. 

Example 6.1. We consider the convection-diffusion (1.1) with = (0,1)^, / = 0, 
w = (1 -I- (a; -I- 1)^/4,0), and the following Dirichlet boundary conditions: 

f u{x, 0) = 1 a; G [0,1], 

[ u{x, 1) = 0 a; G [0, ij. 

With these data, equation (2.4) reduces to eTiU + eUT 2 -I- (4>ii?i){7 = F, that is to 
the following Sylvester equation, 

{eTi + + eUT2 = F, (6.1) 


where the right-hand side F contains the Dirichlet boundary conditions (see section 3). 
We solve this matrix equation using the lyap function in Matlab, from the Control 
Toolbox, and its “vectorized” version (i.e. its Kronecker formulation) by means of 
a sparse direct solver (Matlab backslash operation). The results are reported in 
Table 6.1(left), for a variety of mesh dimensions and viscosity values; the same number 
of grid nodes, and Uy, was used in the x and y directions, respectively. The numbers 
show the superiority of the use of the (dense) Sylvester solver, compared with that of 
a general-purpose sparse direct solver. 

The right portion of Table 6.1 reports the same type of experiments for the 
rectangular domain il = (0,1) x (0,2). A grid consisting of small squares was still used, 
so that a different number of nodes in the two directions was imposed. Both methods 
are sensitive to the unbalanced dimension growth in the two directions, although the 
matrix oriented approach is still faster. Note that for the largest grid, the performance 
of the dense solver lyap deteriorates more significantly. For such large problems, an 
iterative solver should be preferred. 


e 

rix 

Uy 

Vect. 
CPU time 

Matrix 
CPU time 

0.0333 

65 

130 

0.04 

0.025 

0.0333 

129 

258 

0.15 

0.081 

0.0333 

257 

514 

0.76 

0.415 

0.0333 

513 

1026 

4.28 

2.410 

0.0333 

1025 

2050 

30.25 

24.942 

0.0167 

65 

130 

0.03 

0.016 

0.0167 

129 

258 

0.15 

0.093 

0.0167 

257 

514 

0.70 

0.435 

0.0167 

513 

1026 

5.18 

2.376 

0.0167 

1025 

2050 

35.31 

28.691 

0.0083 

65 

130 

0.03 

0.018 

0.0083 

129 

258 

0.16 

0.087 

0.0083 

257 

514 

0.88 

0.631 

0.0083 

513 

1026 

5.15 

3.245 

0.0083 

1025 

2050 

29.99 

24.735 


Table 6.1 


e 

rix 

Vect. 

Matrix 



CPU time 

CPU time 

0.0333 

65 

0.02 

0.008 

0.0333 

129 

0.07 

0.030 

0.0333 

257 

0.37 

0.159 

0.0333 

513 

2.01 

0.898 

0.0333 

1025 

10.95 

6.389 

0.0167 

65 

0.02 

0.005 

0.0167 

129 

0.08 

0.029 

0.0167 

257 

0.39 

0.155 

0.0167 

513 

1.91 

0.899 

0.0167 

1025 

10.95 

6.354 

0.0083 

65 

0.01 

0.005 

0.0083 

129 

0.08 

0.032 

0.0083 

257 

0.38 

0.154 

0.0083 

513 

1.92 

0.891 

0.0083 

1025 

11.03 

6.443 


Example 6.1. Total CPU time for a sparse direct solver and the function lyap, as the viscosity 
and the problem dimension change. Left: square domain. Right: rectangular domain with same 
mesh size and different number of nodes in x and y directions. 


Example 6.2. We consider solving the algebraic problem stemming from the 
equation of Example 6.1 (in fl = (0,1)^) by means of an iterative solver. In the case 
of the linear system (the Kronecker version of (6.1)), both MI20 and AGMG were 

















matrix-equation strategies for convection-diffusion equations 


13 


e 

nx 

FGMRES-I-AGMG 

GMRES-I-MI20 

KPIK 



CPU time (# its) 

CPU time (# its) 

CPU time (# its) 

0.0333 

129 

0.1649 (13) 

0.2843 ( 7) 

0.2131(24) 

0.0333 

257 

0.3874 (15) 

0.5715 ( 8) 

0.2817(32) 

0.0333 

1025 

11.4918 (20) 

8.4540 ( 8) 

1.5002(44) 

0.0333 

1200 

12.5441 (17) 

8.7843 ( 7) 

1.9722(46) 

0.0167 

129 

0.2150 (14) 

0.2750 ( 7) 

0.1638(22) 

0.0167 

257 

0.4356 (14) 

0.6533 ( 9) 

0.3628(32) 

0.0167 

513 

2.0712 (15) 

2.2171 ( 9) 

0.6324(38) 

0.0167 

1025 

11.5428 (18) 

8.0454 ( 8) 

2.8454(64) 

0.0167 

1200 

13.2109 (18) 

9.5501 ( 8) 

2.1961(52) 

0.0083 

129 

0.1501 (14) 

0.2685 (10) 

0.1394(22) 

0.0083 

257 

0.3651 (15) 

0.5871 ( 8) 

0.2885(34) 

0.0083 

513 

1.6615 (14) 

2.1814 (10) 

0.5439(42) 

0.0083 

1025 

10.0859 (18) 

10.6729 (11) 

2.7800(66) 

0.0083 

1200 

14.4866 (18) 

11.0856 ( 9) 

2.8459(58) 


Table 6.2 

Example 6.2. Total CPU time and number of iterations for solving the Sylvester equation by 
iterative methods. 


considered as preconditioners. As discussed in section 4.1, the iterative solver for the 
Sylvester equation (6.1) was the extended Krylov subspace method [27] (KPIK in the 
following). In both solvers, the stopping criterion was based on the relative residual 
2-norm, with stopping tolerance equal to 10“®. Numerical results are reported in 
Table 6.2 as the viscosity and number of grid nodes in each direction change. We 
observe that the number of iterations varies for all methods as the mesh is refined, 
and in a more significant manner for KPIK. However, we recall from section 4.1 that 
memory requirements in this case are not an issue, as the extended Krylov subspace 
actually generated is a subset of R"®. For the two algebraic multigrid preconditioners, 
a mildly varying number of iterations with respect to rix has been largely observed in 
the literature, at least experimentally. All methods are not very sensitive to changes 
in viscosity. In terms of CPU times, the method based on the matrix equation shows 
the best performance, especially as the problem size increases. In summary, iteratively 
solving the problem by resorting to its matrix equation formulation pays off in this 
case. A comparison with the timings of Table 6.1 is also of interest. Note that for 
the largest size, the iterative solvers should be preferred, whereas for most other sizes, 
especially in the linear system formulation, the direct solver is superior. This fact is 
typical of two-dimensional problems, for which sparse direct solvers remain attractive 
over iterative ones, also for very fine discretizations when the Kronecker formulation 
is used. 

In the next two examples we consider problems where the convection coefficients 
both depend on both variables, although in a separable manner. The matrix-equation- 
based strategy is then used as a preconditioner, in which the original operator is 
approximated via a Sylvester operator. 

Example 6.3. We consider the convection-diffusion equation in H = (0,1)^ with 
w = (2?/(l - a;^), -2a;(l - y^)), 

and the Dirichlet boundary conditions u(0,y) = 0, u(l,y) = 0, u(x,0) = 0 and 
u(x, 1) = 1. This is a simple model for the temperature distribution in a cavity with 
a ‘hot’ external wall ({1} x [0,1]) and the wind characterized by w determines a 
recirculating flow; this model was used as a test example in [14]. Data were generated 
with the IFISS package [12], which uses a uniform finite element discretization on 
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a rectangular grid. The coefficient matrices of the Sylvester operator preconditioner 
were generated by finite differences and the corresponding equation solved with KPIK. 
Thanks to the chosen uniform discretization in IFISS, the obtained preconditioning 
operator was still effective, in spite of the different discretization technique. Numerical 
results are reported in Table 6.3, where a stopping tolerance of 10“® was used. The 
preconditioner MI20 broke down several times on this example, therefore its results 
where not included. Instead, the last column shows timings for the Matlab default 
sparse solver (backslash operator). 

The numbers in the table show a very good performance of the new approach 
with respect to the number of iterations, but timings are in general higher than with 
AGMG preconditioning; better performance of FGMRES+KPIK may be obtained 
with a compiled code. Note that both iterative solvers are largely superior to the 
sparse direct method. 


e 

rix 

FGMRES-I-AGMG 

FGMRES-I-KPIK 

DIREGT 



CPU time (# its) 

CPU time its) 

CPU time (# its) 

0.0050 

129 

0.0377 ( 6) 

0.1301 (10) 

0.0785 

0.0050 

257 

0.1216 ( 6) 

0.2509 ( 9) 

0.4338 

0.0050 

513 

0.5970 ( 6) 

0.7709 ( 8) 

3.3971 

0.0050 

1025 

3.1946 ( 7) 

6.0621 ( 8) 

19.0872 

0.0025 

129 

0.0634 ( 7) 

0.1306 (10) 

0.0824 

0.0025 

257 

0.1073 ( 5) 

0.2173 ( 9) 

0.4223 

0.0025 

513 

0.5029 ( 5) 

0.7151 ( 8) 

2.4794 

0.0025 

1025 

2.5261 ( 6) 

5.3601 ( 7) 

16.1251 

0.0013 

129 

0.1882 (10) 

0.1686 (10) 

0.1117 

0.0013 

257 

0.2531 ( 7) 

0.2283 ( 9) 

0.4317 

0.0013 

513 

0.5075 ( 5) 

0.7088 ( 8) 

3.2329 

0.0013 

1025 

2.1558 ( 5) 

3.9864 ( 6) 

24.6895 


Table 6.3 

Example 6.3. CPU time and number of iterations for different preconditioners and a sparse 
direct solver, as the viscosity and mesh parameter vary. 


Example 6.4. We consider a variant of Example V in [11] in Q = (0,1)^, where 
the convection vector was modified so as to obtain a divergence free convection term, 
and so as not to have zero mean, namely 

w = (y(l - (2x + 1)2), -2(2a: + 1)(1 - y^)), 

together with zero Dirichlet boundary conditions except for the side y = 0 where 

f u(a:, 0) = 1 + tanh[10 + 20(2a: — 1)], 0 < x < 0.5, 

( u(x, 0) = 2, 0.5 < X < 1. 

The problem was discretized using centered finite differences, and the grid was selected 
fine enough so as not to have spurious oscillations. As in the previous example, we 
compare the performance of various preconditioners, as the viscosity and the mesh 
parameter change. Results are reported in Table 6.4, where the stopping tolerance 
10 “® was used. 

This is a recognized hard problem due to the strong layer appearing near the 
boundary y = 0. This appears to be the only problem where FGMRES+AGMG gives 
much worse results than the other methods, both in terms of variability on the number 
of iterations, and in terms of CPU time. On the other hand, MI20 and KPIK worked 
very well as preconditioners, with an either constant or even decreasing number of 
iterations, and lower CPU times, with somewhat lower values for the matrix-equation- 
based approach. 
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e 

rix 

FGMRES-I-AGMG 

GMRES-I-MI20 

FGMRES-I-KPIK 



CPU time (# its) 

CPU time [ff its) 

CPU time (# its) 

0.1000 

128 

0.1081 (11) 

0.1989 ( 4) 

0.5743( 8) 

0.1000 

256 

0.3335 (14) 

0.3812 ( 4) 

0.5351( 7) 

0.1000 

512 

1.0773 (11) 

1.1731 ( 4) 

1.0543( 7) 

0.1000 

1024 

9.2493 (17) 

4.3287 ( 4) 

2.6372( 6) 

0.1000 

2048 

52.0430 (15) 

19.6757 ( 4) 

16.7394( 5) 

0.0500 

128 

0.0936 ( 9) 

0.2269 ( 4) 

0.5168(10) 

0.0500 

256 

0.2897 (12) 

0.3862 ( 4) 

0.6455( 9) 

0.0500 

512 

1.2603 (11) 

1.2380 ( 4) 

1.1769( 8) 

0.0500 

1024 

10.3623 (18) 

4.3345 ( 4) 

3.0812( 7) 

0.0500 

2048 

60.5041 (17) 

20.4056 ( 4) 

14.9237( 6) 

0.0333 

128 

0.0882 ( 8) 

0.2368 ( 4) 

0.6428(11) 

0.0333 

256 

0.2181 ( 9) 

0.4218 ( 4) 

0.8149(11) 

0.0333 

512 

1.5849 (14) 

1.1977 ( 4) 

1.3786( 9) 

0.0333 

1024 

5.4624 (12) 

4.4130 ( 4) 

3.7214( 8) 

0.0333 

2048 

120.9686 (23) 

20.1120 ( 4) 

17.9188( 7) 


Table 6.4 

Example 6.4- CPU time and number of iterations for different preconditioners, as the viscosity 
and mesh parameter vary. 


6.2. Three-dimensional problems. Like in the previous section, we first re¬ 
port on the case where the discretized problem can be directly solved as a Sylvester 
equation, and then on the case when the Sylvester equation serves as preconditioner 
for the more involved problem. 

Example 6.5. We consider the 3D problem —eAu -I- w • Vu = 1, in ft = (0,1)^, 
with convection term 


w = (xSinx, y cosy, e ), 

and zero Dirichlet boundary conditions. After a centered finite difference discretiza¬ 
tion with Ux = riy = Uz nodes in each direction, the equation (5.2) takes the form 

{I(g>Ti)U+U T3 + {t4 ®/) U+{I(g>^iBi)U+l4 B3T3 + [{'$2B2f <S>I]U = 11^, 

where 1 is the vector of all ones, and this corresponds to the following Sylvester 
equation in U, 

[I ® (Ti + $iBi) + {T 2 + ■^ 2 B 2 f ®I]U+U{T3 + B 3 T 3 ) = 11^. (6.2) 

Note that the fact that the forcing term in the original equation yields a low rank 
right-hand side in the matrix equation is crucial for the solution process. 

We solved this matrix equation by means of KPIK, and compared its performance 
with the solution of the corresponding linear system (via its Kronecker formulation) 
with GMRES preconditioned by MI20, and with flexible GMRES preconditioned by 
AGMG. Results are reported in Table 6.5, with stopping tolerance equal to 10“®. 

As in the corresponding 2D problem (Table 6.2), the number of iterations is quite 
stable for both algebraic multigrid preconditioners, whereas it varies, though very 
mildly, when KPIK is used. On the other hand, CPU time is very much in favor of 
the matrix equation solver, whose computation exploited the lower order complexity 
of the problem. One order of magnitude lower timings can be observed in some 
instances. All strategies seem not to be sensitive to the choice of viscosity in the 
considered range. 
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e 

rix 

FGMRES-I-AGMG 

GMRES-I-MI20 

KPIK 



CPU time its) 

CPU time (7^ its) 

CPU time (# its) 

0.0050 

50 

0.6044 (13) 

1.0591 ( 6) 

0.1353 (18) 

0.0050 

60 

1.1250 (14) 

1.8022 ( 7) 

0.1734 (18) 

0.0050 

70 

2.0385 (14) 

3.4253 ( 7) 

0.2326 (20) 

0.0050 

80 

3.4803 (14) 

4.4297 ( 7) 

0.3583 (20) 

0.0050 

90 

5.7324 (15) 

6.8705 ( 7) 

0.4999 (22) 

0.0050 

100 

8.0207 (15) 

9.7207 ( 7) 

0.5677 (22) 

0.0010 

50 

0.6556 (13) 

1.1306 ( 6) 

0.1811 (18) 

0.0010 

60 

1.3011 (14) 

1.7854 ( 7) 

0.2386 (18) 

0.0010 

70 

1.9509 (14) 

2.7829 ( 7) 

0.2346 (20) 

0.0010 

80 

3.5291 (14) 

4.6576 ( 7) 

0.4096 (20) 

0.0010 

90 

5.1344 (14) 

6.8176 ( 7) 

0.4253 (22) 

0.0010 

100 

7.6815 (14) 

9.4935 ( 7) 

0.5446 (22) 

0.0005 

50 

0.7039 (14) 

1.0530 ( 6) 

0.1751 (16) 

0.0005 

60 

1.2560 (14) 

1.7341 ( 6) 

0.2314 (18) 

0.0005 

70 

2.2242 (14) 

2.9667 ( 7) 

0.2301 (20) 

0.0005 

80 

3.4558 (14) 

4.5964 ( 7) 

0.3472 (22) 

0.0005 

90 

4.8076 (14) 

6.4841 ( 7) 

0.4257 (22) 

0.0005 

100 

7.3914 (14) 

9.6274 ( 7) 

0.5927 (24) 


Table 6.5 

Example 6.5. CPU time and number of iterations for different preconditioners, as the viscosity 
and mesh parameter vary. 


Example 6 .6. Finally, we consider the 3D convection-diffusion equation in = 
(0,1)^, with / = 1 and convection vector 

w = (yz(l-x2),0,e^), 

and zero Dirichlet boundary conditions. Due to the structure of w, we can treat 
separately the (z) variable and the {x,y) variables^. With this variable ordering in 
mind, we can define U of size x (nxUy), and the following corresponding matrix 
equation 

TsU+U{I^T2+tI + + = 11^. (6.3) 

With the approximation Ti « vil, the matrix-equation-based preconditioner is given 

by 

P : V ^ {T 3 + T 3 B 3 )V + V{I( 8 >T 2 + T^ (E)I + I<S> vi^iBi), (6.4) 

with coefficient matrices of dimension x riz and UxUy x UxTiy. We compare the 
performance of flexible GMRES with this preconditioner (the Sylvester equation is 
iteratively solved as described in section 4.1), with that of flexible GMRES precon¬ 
ditioned by AGMG, and of GMRES with MI20, both applied to the corresponding 
linear system. The stopping tolerance was set to 10“®. The results of our numeri¬ 
cal experiments are shown in Table 6.6. We readily notice that GMRES-I-MI20 did 
not perform well, and we had to stop the solution process for the two smaller values 
of the viscosity, due to the excessive GPU time necessary for building the precondi¬ 
tioner. The other preconditioning strategies perform comparably all the way up to 10® 
unknowns, with a steadily lower number of iterations for the matrix-equation-based 
preconditioner, which in turns implies lower memory requirements of flexible GMRES. 

^ Other variable aggregations are possible. The one we chose allowed us to explicitly treat the 
first derivative in the z direction in the preconditioner. 
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In terms of CPU time performance is comparable, although the new strategy usually 
requires lower time; we expect that the use of a compiled implementation of this new 
strategy will make this difference more remarkable. 


e 

fix 

FGMRES-I-AGMG 

GMRES+MI20 

FGMRES-I-KPIK 

0.5000 

50 

0.6239 (15) 

0.9182 ( 6) 

0.9291 ( 6) 

0.5000 

60 

1.2095 (15) 

1.8236 ( 7) 

1.2027 ( 6) 

0.5000 

70 

2.1041 (15) 

3.1649 ( 7) 

1.6585 ( 6) 

0.5000 

80 

3.7370 (16) 

4.9765 ( 7) 

2.4943 ( 6) 

0.5000 

90 

7.5874 (16) 

9.2040 ( 8) 

3.2513 ( 6) 

0.5000 

100 

7.7626 (16) 

11.9912 ( 8) 

4.7548 ( 6) 

0.1000 

50 

0.8041 (17) 

75.8874 ( 6) 

1.4610 ( 8) 

0.1000 

60 

2.1310 (18) 

- (-) 

1.5299 ( 8) 

0.1000 

70 

2.8043 (18) 

- (-) 

1.8926 ( 8) 

0.1000 

80 

5.1219 (19) 

- (-) 

3.2928 ( 9) 

0.1000 

90 

7.3179 (19) 

- (-) 

4.6429 ( 9) 

0.1000 

100 

9.5759 (19) 

- (-) 

6.5590 ( 9) 

0.0500 

50 

0.6780 (17) 

158.5447 ( 6) 

1.1997 (10) 

0.0500 

60 

1.4318 (18) 

- (-) 

1.7296 (10) 

0.0500 

70 

2.8427 (19) 

- (-) 

2.5215 (10) 

0.0500 

80 

4.9616 (20) 

- (-) 

3.6615 (10) 

0.0500 

90 

7.1038 (20) 

- (-) 

5.0098 (10) 

0.0500 

100 

10.7181 (21) 

- (-) 

6.8661 (10) 


Table 6.6 

Example 6.6. CPU time and number of iterations for different preconditioners, as the viscosity 
and mesh parameter vary. stands for excessive CPU time in building the preconditioner. 


7. Conclusions and outlook. We have implemented a new matrix-equation- 
based strategy for solving or preconditioning a variety of two and three-dimensional 
convection-diffusion problems. Our preliminary numerical experiments show that 
the new approach performs comparably well with respect to state-of-the-art alge¬ 
braic multigrid preconditioners. As opposed to earlier attempts with this type of 
approaches, we have shown that recently developed Sylvester equation solvers can 
make the matrix equation strategy very appealing whenever the original partial dif¬ 
ferential equation has separable coefficients. 

Our current implementation is limited to the use of uniform meshes for rectangu¬ 
lar or parallelepipedal domains. We plan to generalize the approach to more general 
settings in the future. The restriction to separable convection terms can be relaxed 
by appropriately approximating the operator at the preconditioning level, by design¬ 
ing special corresponding separable approximations. These possibilities will also be 
explored. 

Our description could be generalized to tensor structures with more than three 
dimensions, as they occur in many emerging applications [16],[17],[18]. We will explore 
the possibility of extending our algorithmic methodology to this setting, taking into 
account the role of the convection terms in this more advanced framework. 
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